Towards
Optimizing Explicit Fracture Simulations of Unconventional Reservoirs
Daniel
Barreca
Summary
Simulations
are carried out using MRST’s black oil module on reservoirs with properties
similar to the Bakken Shale. Simulations of a fractured reservoir based on
fracture simulations are also carried out. Hydraulic fractures are treated
through the explicit fracture modelling method. In an attempt to reduce grid
complexity, the effect of fracture geometry and fracture interaction is first
investigated. Results demonstrated that production is dependent on fracture
geometry, therefore the mesh cannot be simplified to structured meshed for
explicit fracture systems. Additionally,
the method of snapshots is applied to one of the reservoirs. Results from computing
the basis functions and eigenvalues from this
snapshot matrix demonstrate the potential of reduced order modelling on
explicitly fractured unconventional simulations.
Introduction
Due to the large amount of oil and gas held
in unconventional reservoirs, it is important that efficient and accurate
models exist for simulation of these resources. When trying to model
unconventional reservoirs, certain challenges arise such as matrix/fracture
interactions, nonlinear pressure distributions, fracture networks, and geomechanical considerations (Farah 2016).
Modelling complex physical processes such as these may be computationally
demanding, facilitating a need for methods that reduce CPU time of simulation.
One way of reducing complexity of a problem is through a reduced order model
(ROM).
Reduced order modelling is
used to reduce computational demands of a simulation by reducing the degrees of
freedom of a problem. Various ROM techniques exist, though the focus of this
project is using a classification of ROMs known as Proper Orthogonal
Decomposition (POD). More specifically, to attempt to reduce the computational
demands of unconventional reservoir simulation, the goal of this project is to
test the applicability of the method of snapshots to such reservoirs. The
method of snapshots, developed by Sirovich, is a
technique of producing a POD basis for a ROM (Sirovich
1987). POD methods have been used to
create ROMs for simple reservoir simulations in a number of studies (Jing et
al. 2017, Cardoso et al. 2009, Markovinovic et al.
2002). Other POD studies more specifically related to unconventional simulation
include ROMs for geomechanics and coupled flow (Florez
and Gildin 2018), multiscale fracture simulation (Oliver
et al. 2017), and thermal simulation of oil shale (Rousset
2010). For this project, explicitly fractured simulations are first simulated
and compared. The method of snapshots will then be applied explicitly fractured
systems, as well as be applicable to simulation of my chosen reservoir
(Bakken).
Hydraulic
Fracturing
Reservoirs such as shale have extremely low
values of porosity and permeability. Hydraulic fracturing is a technique where
fluid is injected into a rock formation and proppants are embedded, allowing
for a highly permeable flow path towards the wellbore. This method has made
unconventional resources an attractive source of economic value. As a result,
simulation of these resources (along with hydraulic fractures), must be
accurately calculated. As one of the steps in building the reservoir model, hydraulic
fracture pathway must be considered since this fracture will need to be
modelled.
Though many models assume
the fracture will extend linearly, the hydraulic fracture geometry may not
always extend perpendicular and linearly from the wellbore. Various factors can
generate more complex fracture geometries. For example, in cases where the
reservoir is naturally fractured, the hydraulic fracture may follow along this
pathway. Another consideration for fracture geometry is the stress field within
the reservoir. There is a possibility of hydraulic fractures interacting with
each other due to created stresses and depending on the fracture spacing (Wu
2014). Additionally, fracture path may be influenced by the pre-existing stress
field, as indicated in fracture simulations by Wu (2014), shown in Figure 1.
Hydraulic fracturing possibilities such as these will be investigated through
reservoir simulation.

![]()
Preprocessing
Procedure for Simulating More Complex Geometries
Because of the geometric complexity of
hydraulic fractures, more advanced meshing and geometry generation software
packages are used for preprocessing. Autodesk INVENTOR is a commercial 3D CAD
software that is commonly used in various industries for CAD creation. The
software package is used here due to its flexibility and advanced geometry
creation features. Various fracture models can be created, including complex
fracture geometries. A simple example is shown in Figure 2.
The CAD
file is exported as a (.dwg) file and imported into ANSYS ICEM CFD. ICEM CFD is
a commercial mesh generation software package. Again, a different software is
used here for meshing because ICEM CFD is more flexible than MATLAB in
importing CAD files. In addition to this, ICEM CFD has more advanced meshing
control and generation options than MRST. Cell sizings,
cell growth rates, and various meshing algorithms can easily be implemented to
mesh the geometry. An example of a generated triangular mesh is again shown in
Figure 2. Cell sizing can be straightforwardly refined or coarsened by
specifying a cell size or number of nodes around specified curves, allowing for
controlled node density and growth rates around hydraulic fractures.
From
here, the mesh can be exported into specific formats for various simulation
packages such as FLOW3D or ANSYS FLUENT. MRST is not one of the available
options for mesh exporting. As a workaround, the mesh is exported as a (.msh) text file, which is the unstructured mesh format that
ANSYS FLUENT accepts. X and Y coordinates from the nodal positions are
extracted from this text file and stored as vectors in MATLAB. Finally, using
these two vectors, the mesh can be reconstructed using MRST’s TriangularGrid() function.
Fracture permeability and porosity cannot simply be
assigned the correct values based on location because fracture geometry may be
complex. Instead, to assign permeability and porosity to the fracture, all
cells are looped over and the volumes are checked. A histogram can be created
to determine which size cells lie within the fracture. Since the mesh is most
refined in these sections, cell volumes should be smallest within the
fractures. Unfortunately, this is not always the case. As shown in Figure 2,
using this method will result in slight deviations from the true fracture, as
random cells outside, but near, the fracture may be included. However, a decent
representation of the fracture using this method is attained. A more accurate
method of assigning correct permeability values would include some
implementation where cell centroid locations are utilized.

Fig.
2
Preprocessing Procedure
General
Simulation and Geometry Details
Because the Bakken shale is an oil shale, the
simulations are carried out using MRST’s black oil module. All simulations are
completed using a 2D discretization due to the computational requirements of
the 3D meshes (further explained in a later section). The geometries for all
following simulations have the dimensions of one of the figures specified in
Figure 4. For the simulations with 3 hydraulic fractures, a fracture spacing of
250 ft is used, which is based off common spacings from fracturing done in the
Bakken Shale. For the simulations consisting of 4 hydraulic fractures, a
spacing of 150 ft is sued, which is based off simulations from Wu (2014). Most
reservoir and fluid properties are based off the Bakken Shale and specified in
the literature review. Properties not clearly defined or not included in the
literature review are included in Figure 3.
|
Property |
Value |
Unit |
|
Matrix Porosity |
6 |
% |
|
Matrix Permeability |
0.01 |
mD |
|
Fracture Porosity |
20 |
% |
|
Fracture Permeability |
10 |
mD |
|
Fracture Width |
0.2 |
in |
|
Wellbore Radius |
0.4 |
in |
|
Wellbore Pressure |
1000 |
psi |
|
Simulation Time |
10 |
years |
|
Simulation Timestep |
0.0025 |
years |
Fig.
3
Extra Details for
Simulation, not included in Literature Review


Fig.
4
Base Geometry Dimensions
Grid
Independence Study
To ensure special convergence of the
simulations, a grid independence study is initially carried out. Because
pressure drawdowns, fluid properties, and rock properties are similar between
all geometries, only one grid independence study is completed. The grid
independence study is run this single geometry, and a similar global mesh sizing
is to be used for the rest of the geometries. Three meshes of different node
densities are created for the straight-line perforations’ geometry. 1000, 500,
and 250-feet cell-sizings are used for the global
mesh size, resulting in 14000, 7000, and 5000 nodes (Fig. 5). The primary value
of interest, oil produced, is monitored and recorded over a 10 year simulation
period. The results are displayed in Figure 6. While a significant difference
in production exists between the 5000 node and the 7000 node system, this
difference is insignificant between the 7000 and 14000 node systems. As a
result, a global grid spacing of 500 ft will be specified in all following
simulations.

Fig.
5
Grid Independence
Study: 5000 nodes (left), 7000 nodes (middle), and 14000 nodes (right)

Fig.
6
Grid Independence
Study, Monitoring Cumulative Oil Production
Impact
of Fracture Geometry
Before moving on, an initial and important
question is whether fracture geometry will have a significant impact on
production. If fracture geometry does not have an impact on production, then
there is little point in modelling anything beyond simple, straight hydraulic
fractures. Straight hydraulic fractures have an inherent advantage for
simulation in that they can be modelled easily on structured meshes. Structured
meshes have various advantages such as computational efficiency and alignment
with flow or permeability fields. Setting up a structured mesh for a more
complex hydraulic fracture system (such as that in Figure 2) would be at the
least very difficult to set up and a manually labor-intensive process. In many
cases, structured meshes may even be practically impossible to generate.
Two
situations are tested in this study: effect of fracture curvature (Fig. 7) and
effect of fracture interaction (Fig. 8). For the fractures in Figure 7, it was
ensured during the geometry creation process that every fracture had a
half-length of 400 ft, spacing of 250 ft, and node density throughout the mesh.
Thus, all fractures have an equivalent pore volume, and the only difference is
the degree of curvature. For the fractures in Figure 8, the effect of fracture
interaction is simulated by roughly replicating a geometry from Wu (2014).
Simulations are run for 10 years at a
timestep of 0.0025 years with a wellbore pressure of 1000 psi. Cumulative
production is plotted for each situation in figures 9 and 10. Additionally, pressure
contour plots are generated and displayed in figures 11 and 12. The results
demonstrate that fracture geometry may have an impact on cumulative production.
In the case of analyzing
fracture curvature, a significant curvature decreases production, while a
slight curvature increases production. The contour plots suggest that an
increase in curvature results in an increase in near wellbore production. This
can be seen at a timestep of 1 year, where pressure in areas surrounding the
wellbore is lower for higher curvature fractures. This suggests that there may
be some optimal curvature or angle that optimizes production. Further analysis
would need to be completed to fully investigate this possibility.
In the case of analyzing
fracture interaction, the results are as expected. When considering fracture
interaction, the system with shortened inner fractures results in a noticeably
lower production value. As demonstrated in the contour plots, this is simply
due to longer fractures permeating further into the reservoir, and consequently
draining more of the reservoir. One factor that was not considered here was
that inner fractures merging with the main fractures may extend the main
fractures’ length. Further literature review needs to be completed on this
topic to fully understand hydraulic fracture interaction.

Fig.
7:
Varying Hydraulic
Fracture Curvature

Fig.
8
Varying Hydraulic
Fracture Interaction

Fig. 9
Cumulative Production Plot of Results from Simulations
in Fig. 7

Fig. 10
Cumulative Production Plot of Results from Simulations
in Fig. 8

Fig.
11
Pressure (psi)
Contour Plots for 5 Time Snapshots: Effect of Curvature

Fig.
12
Pressure (psi)
Contour Plots for 5 Time Snapshots: Effect of Fracture Interaction
Reduced
Order Modelling in Reservoir Simulation
While simpler reservoir and fracture
simulations such as these may not require a large amount of computational
resources, situations such as 3D simulations, realistic sized reservoirs, and
increased number of fractures may take a significant amount of time to
complete. Take Figure 13 as an example. When no fracture is modelled, the mesh
only consists of 710 nodes. After a fracture is added in and discretized with
the reservoir, the node count drastically increases to 44,913 nodes, which is
63 times the number of nodes in the standard simulation. This also does not
take into consideration the mesh shown is a 2D system. When extending to 3D, it
may not be as simply as creating a few layers. Since the cell sizes around the fractures are on the order of magnitude of
inches, the vertical resolution should also be this order of magnitude in order
to maintain favorable aspect ratios and other important cell quality criteria.

Fig.
13
Discretization
Increase When Including Fracture
In addition to a drastic
increase in cell count, an engineer may want to run hundreds or thousands of
simulations related to production optimization or sensitivity analysis. The
question proposed is how this high dimensional system can be solved in a more
efficient manner. A potential solution to this question is reduced order
modelling. There are various methods of reduced order modelling, one of those
being Proper Orthogonal Decomposition (POD) based methods.
The time evolution of pressure
depends on a nonlinear function of variables such as space, time, and rock
properties.
![]()
Normally, using MRST, a finite volume discretization
of this system will be carried out, which is a finite representation of the
infinite dimensional space. Essentially, the PDE is substituted for a large
system of solvable differential equations. Because explicit fracture systems
significantly increase the discretization of the reservoir system, we would
like to reduce this large system in some way. The first step in the reduced
order modelling process is to assume that the spatial (x) dependence
and the time (t) dependence of a problem can be separated through the
separation of variables:
![]()
The functions of interest that need to be
determined in the equation above are the basis functions,
.
Afterwards, a finite sum of basis
functions may be determined based on the time dynamics of the system.
To
determine these basis functions, the full MRST simulation must be run for a
number of time steps so that a snapshot matrix for the variables of interest
can be attainted. For example, the snapshot matrix, P, is simply just a
collection of the pressure values of each grid cell at time interval ‘NT’. This
information is stored in an (n x NT) matrix where n is the number of grid
cells.
![]()
After this matrix is stored, an SVD may
be computed on the matrix set.
![]()
The matrix U, which is computed as part of
SVD, contains information about the dominant modes for the basis
functions.
The
process described above is completed for the fracture interaction system from
Figure 13. Eigenvalues are plotted in Figure 14 and basis function energy
contour plots are shown in Figure 15. This system may be suitable for reduced
order modelling, as indicated by Figure 14. The eigenvalues in figure 14
demonstrated that the energy of the system varies by
orders of magnitudes. Therefore, the vast
majority of energy in this system can be represented by a small number of modes.
The energy of these modes can be viewed and compared, shown in Figure 15. The
first 8 basis functions(which are the most dominant energy systems) are vastly
different than, for example, the 200th largest energy mode (which is
essentially just noise).

Fig.
14
Semilog plot of Eigenvalues for the
Fracture Interaction System

Fig.
15
Basis Functions of
the Complex Fracture Simulation
Conclusions
1.)
Modelling accurate fracture geometries during
the simulation process is important for accurate estimation of resources, Fracture
geometry may have an impact on production. In some cases, fracture curvature
can increase production, though a large curvature may significantly decrease
production. Hydraulic fracture interaction should also be considered.
2.)
Modelling fractures through an explicit
method will significantly increase cell count. Reduced order modelling of
unconventional has potential, as indicated by an SVD on the system.
Future
Work
1.)
More work needs to be done in determining
effect of fracture geometry. A larger number of fractures, different fracture
spacings, various reservoir properties or fluid properties, and fracture
lengths in comparison to reservoir length should be investigated to fully
understand the impact of fracture geometry.
2.)
Different scenarios of fracture interaction
may be considered to estimate the impact on production
3.)
Reduced order modelling (either non-intrusive
or intrusive) should be carried out. The results in this study simply show the
potential of reduced order modelling, but a working reduced order model was
never complete due to possible errors in the code or lack of understanding of
the problem.
References
|
Farah, N., 2016. Flow Modelling in Low Permeability
Unconventional Reservoirs (Doctoral dissertation, Paris 6). |
|
Sirovich, L., 1987. Turbulence and the dynamics
of coherent structures. I. Coherent structures. Quarterly of applied
mathematics, 45(3), pp.561-571. |
|
Cao Jing, Zhao Hui,
Yu Gaoming, and Xie Yunhong, 2017. The POD Model Order Reduction Method
of Accelerating Reservoir Numerical Simulation. IOSR Journal of Applied Geology and Geophysics, Volume 5, Issue 1 Ver. I, pp. 107-114. |
|
Cardoso, M.A., Durlofsky,
L.J. and Sarma, P., 2009. Development and
application of reduced‐order modeling procedures for
subsurface flow simulation. International journal for numerical
methods in engineering, 77(9), pp.1322-1350. |
|
Markovinovic, R., Geurtsen,
E.L., Heijn, T. and Jansen, J.D., 2002, September.
Generation of Low-Order Reservoir Models Using POD, Emperical
Grammians and Subspace Identification. In ECMOR
VIII-8th European Conference on the Mathematics of Oil Recovery. |
|
Florez, H. and Gildin,
E., 2018, September. Model-Order Reduction Applied To Coupled Flow and
Geomechanics. In ECMOR XVI-16th European Conference on the
Mathematics of Oil Recovery. |
|
Oliver, J., Caicedo,
M., Huespe, A.E., Hernández, J.A. and Roubin, E., 2017. Reduced order modeling strategies for
computational multiscale fracture. Computer Methods in Applied
Mechanics and Engineering, 313, pp.560-595. |
|
Rousset, M., 2010. Reduced-order
modeling for thermal simulation (Doctoral dissertation, Stanford
University). |
|
Wu, K., 2014. Numerical
modeling of complex hydraulic fracture development in unconventional
reservoirs (Doctoral dissertation). |